M1.D3: Schrodinger’s hydrogen atom

Learning objectives

• When given the time-independent Schrodinger equation identify the different parts of the Hamiltonian. Identify the eigenfunctions and eigenvalues as the solutions of the Schrodinger equation.

• Describe how the quantum numbers and their rules are necessary

• Plot with matplotlib the energy levels of the hydrogen atom and identify their degeneracy

• Plot the radial distribution function of the electron in hydrogen. Explain how we can use it to tell how far the electron is.

• Distinguish between most probable value and expectation (average) value of the electron.

• Plot 3D wavefunctions and coming up with an isosurface

Concepts that need to be learned (remembered?) for the hydrogen atom.

• principal quantum number: \(n=1,2,3,4...\)

• angular or azimuthal quantum number : \(\ell=0,1,2..n-1\)

• magnetic quantum number \(m_l = -\ell,-\ell+1,...,-1,0,1,\ell-1,\ell\)

The Hamiltonian for the Hydrogen Atom

The first step to figure out a quantum system is to solve the time-independent Schrodinger Equation

\[ \hat{H}\Psi = E\Psi \]

We are not going to solve step by step this equation for the Hydrogen atom (check the references for the mathematical development), rather, we are looking at the main points so that students can achieve the learning objectives.

What changes for each system is the Hamiltonian \(\hat{H}\). The Hamiltonian operator has the form

\[\hat{H} = -\frac{\hbar^2}{2m}\frac{\delta^2}{\delta x^2} + V(x)\]

Unlike the particle in a box, now in the hydrogen atom there is an electrostatic attractive potential between the electron and the proton. Considering their charges and the distance between them. $\( V(r) = \frac{1}{4\pi\epsilon_o}\frac{q_{e^{-}} q_{p^{+}}}{r}\)\( The above equation is in metric units, but if we define an energy unit called Hartree (not metric) we can make \)4\pi\epsilon_o = 1\( and then we have \)\frac{q_{e^{-}} q_{p^{+}}}{r}$

The Hamiltonian is $\(\hat{H} = -\frac{\hbar^2}{2\mu}\frac{\delta^2}{\delta (x,y,z)^2} + \frac{q_{e^{-}} q_{p^{+}}}{r} \)$

where \(\mu\) is reduced mass for the proton-electron system \(\mu = \frac{m_e \cdot m_p}{m_e + m_p}\). So the Schrodinger equation is

\[-\frac{\hbar^2}{2\mu}\frac{\delta^2\Psi}{\delta (x,y,z)^2} + \frac{q_{e^{-}} q_{p^{+}}}{r}\Psi = E\Psi \]

Commonly, the first derivative with respect to a multidimensional variable is a vector and it is indicated with the symbol \(\nabla\). So for the first derivative of any function we write what is called the gradient

\[ \nabla f = \frac{\delta f}{\delta x} + \frac{\delta f}{\delta y} + \frac{\delta f}{\delta z} \]

and for the second derivative we actually have more than 3 elements because of all cross-derivatives

\[ \nabla^2 f = \frac{\delta^2 f}{\delta x^2} + \frac{\delta^2 f}{\delta y \delta x} + \frac{\delta^2 f}{\delta z \delta x} + \frac{\delta^2 f}{\delta y^2}... \]

Using this notation the Schrodinger equation is

\[ -\frac{\hbar^2}{2\mu}\nabla^2 \Psi + \frac{q_{e^{-}} q_{p^{+}}}{r}\Psi = E\Psi\]

We won’t be spending any time solving that differential equation, but there are few aspects that need to be kept in mind.

While the potential depends only on the distance between particles “r”, the derivatives are over the three coordinates of the space (x,y,z). This creates a coordinate problem that forces us to convert the x,y,z into three other coordinates called spherical coordiantes(\(r,\theta,\phi\))

The Schrodinger equation in these coordinates may look more complicated

\[ -\frac{\hbar^2}{2\mu}\Big [\nabla^2_r +\frac{1}{r^2}\nabla_{\theta,\phi}^2 \Big ]\Psi(r,\theta,\phi) -V(r)\Psi(r,\theta,\phi)=E\Psi(r,\theta,\phi) \]

but it can actually be solved more easily by splitting it into a radial part (\(R(r)\)) and an angular part (\(Y(\theta,\phi)\)) such that

\[ \Psi(r,\theta,\phi) =R(r)Y(\theta,\phi) \]

So we can rewrite the overall Schrodinger equation for the Hydrogen atom as:

\[ Y(\theta,\phi)\Big [-\frac{\hbar^2}{2\mu}\nabla^2_r -V(r)-E\Big ]R(r) - R(r) \frac{\hbar^2}{2\mu r^2}\nabla_{\theta,\phi}^2 Y(\theta,\phi) = 0 \]

Solutions to the Schrodinger equation for hydrogen: the electronic orbitals

The separation of coordinates in our previous equation will lead into a wavefunction with a radial and an angular solution.

\[ \Psi(r,\theta,\phi) =R(r)Y(\theta,\phi) \]

To translate what is above into English, we mean that the hydrogen wavefunctions (orbitals) will have two parts, a radial part that dictates the size of the orbital (1s,2s,3s…) or how far the electron is from the nucleus (r), and an angular part that dictates the shape of the orbital (s,p,d,f…)

\[ Orbital = RadialSize(r) \times AngularShape(\theta,\phi) \]

• The wavefunction for the radial part \(R_{nl}(r)\) contained what is called the Laguerre polynomials \(L^{2l+1}_{n-l-1}\)

\[ R_{nl}(r) = \sqrt{\Big(\frac{2}{n a_0}\Big)^3 \frac{(n-l-1)!}{2n (n+l)!}} e^{-r/n a_0} \Big( \frac{2r}{na_0}\Big)^l \cdot L^{2l+1}_{n-l-1} \Big(\frac{2r}{n a_0} \Big)\]

• The wavefunction for the angular part \(Y_{lm}(\phi,\theta)\) are called the Spherical Harmonics and they contain the also famous Legendre Polynomials \(P_{lm}(cos \phi)\)

\[ Y_{lm}(\phi,\theta) = \sqrt{\frac{2l+1}{4\pi} \frac{(l-m)!}{(l+m)!} } P_{lm}(cos \phi) \cdot e^{im\theta} \]

the energy of the orbitals (the eigenvalues) unsurprisingly looks like Bohr’s hydrogen atom, just that this time we obtained it by using the Schrodinger equation.

\[ E_n = -\frac{\mu e^4}{8 \epsilon_0^2 h^2 n^2} \]

You can see here the actual functions depending on the values of n l and m

http://chemdata.umr.umn.edu/chem2331/orbitals/h.html

http://chemdata.r.umn.edu/chem2331/rdf/

• Looking at the radial function what values of n and l will write a negative value inside the square root and therefore invalid?

• Looking at the angular function, what values of l and m will write a negative value inside the square root and therefore invalid?

• Write the overall rules of n, l, and m based on your previous answers.

The energy levels in the hydrogen atom

Given the energy from the Schrodinger equation $\( E_n = -\frac{\mu e^4}{8 \epsilon_0^2 h^2 n^2} \)$ If we use SI units for all those magnitudes

• \(\mu = \frac{m_e \cdot m_p}{m_e + m_p}\approx m_e = 9.11x10^{-31} kg\)

• \(e^- = 1.6x10^{-19} C\)

• \(\epsilon_0 = 8.85x10^{-12} F m^{-1}\)

• \(h = 6.63x10{-34} J s\)

\[ E_n = -\frac{2.17x10^{-18} J}{n^2} \; or \; E_n = -\frac{13.6 eV}{n^2} \]

What we also know is the degeneracy

#Lets plot the energy levels applying the quantum number rules
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(15,10))
ax = plt.subplot()

for n in range(1,5):
    ene = -13.6/n**2
    delta = 0
    for l in range(n):
        for m in range(-l,l+1):
            #print("n=",n," l=",l,"ml",m)
            ax.hlines(ene,xmin=delta,xmax=delta+0.5)
            delta = delta + 1
ax.set_title("Energy of Hydrogen orbitals including degeneracy")
ax.set_ylabel("Energy (eV)");

• Look at the code above. What is the value “delta” used for?

• What is the degeneracy level for n=1, n=2, n=3, and n=4? Can you come up with a rule that gives us the degeneracy at a given level “n”?

The radial solution: the size of the atom

If we go back to the the radial part of the wavefunction (\(R_{nl}(r)\)), we see that the Laguerre polynomials \(L^{2l+1}_{n-l-1}(\frac{2r}{n a_0})\) are what makes the wavefunction different at different levels. These polynomials are so well known and ubiquitous that are included in the scipy python package.

\[ R_{nl}(r) = \sqrt{\Big(\frac{2}{n a_0}\Big)^3 \frac{(n-l-1)!}{2n (n+l)!}} e^{-r/n a_0} \Big( \frac{2r}{na_0}\Big)^l \cdot L^{2l+1}_{n-l-1} \Big(\frac{2r}{n a_0} \Big)\]

We will use scipy to plot the radial wavefunction of hydrogen. We will use atomic units: \(a_0=1,\hbar=1,m_{e^-}=1,q_{e^-}=1\)

#Lets plot the radial function
import numpy as np
import scipy.special as spe
def psi_R(r,n,l):
    coeff = np.sqrt((2.0/n)**3 * spe.factorial(n-l-1) /(2.0*n*spe.factorial(n+l)))
    laguerre = spe.assoc_laguerre(2.0*r/n,n-l-1,2*l+1)
    return coeff * np.exp(-r/n) * (2.0*r/n)**l * laguerre

ax1 = plt.subplot(111) 
r = np.linspace(0,10,100)
n=2
l=1 #2p
R = psi_R(r,n,l)
ax1.plot(r,R,lw=2)
ax1.set_title("Radial wavefunction for n="+str(n)+" and l="+str(l))
ax1.set_xlabel("Distance to nucleus in bohrs $a_{0}$");
ax1.set_ylabel("Radial wavefunction $R_{nl}(r)$");

• Spend some time understanding every line in the code above. What does the function psi_R do?

• What variables decide the orbital that we want to plot?

• What is the variable “r” and what is linspace for?

fig = plt.figure(figsize=(12,7))
ax1 = plt.subplot(121) #121 means: the subplots will be 1 row 2 columns and this one is subplot number 1
ax2 = plt.subplot(122, sharey = ax1) #share-y means that it will share the y-axis with the ax1 plot

r = np.linspace(0,10,100)
n=3
l=2 #3d
R = psi_R(r,n,l)
ax1.plot(r,R,lw=2)
ax1.set_title("Radial wavefunction for n="+str(n)+" and l="+str(l))

n=3
l=1 #3p
R = psi_R(r,n=3,l=1)
ax2.plot(r,R,lw=2);
ax2.set_title("Radial wavefunction for n="+str(n)+" and l="+str(l));

Calculate the size of atoms

Remember that the wavefunction plot above does not carry any physical meaning. The square of the wavefunction however gives the probability or probabiilty density.

\[ \text{Probability density}(r,\theta,\phi) = |\Psi(r)\cdot\Psi(\theta,\phi)|^2 \]

Because we are interested in knowing how far from the nucleus the electron is, we could just look at the square of the radial wavefunction.

\[ \text{Radial probability density}(r)=|\Psi(r)|^2 \]

Another magnitude more useful other than the probability is what’s called the radial distribution function where we just consider the probability as a function of r and integrate over all values of \(\theta\) and \(\phi\).

\[ \int_0^\infty\int_0^\pi \int_0^{2\pi} |\Psi(r)|^2 r^2 sin(\theta)dr d\theta d\phi \]

So, integrating for \(\theta\) and \(\phi\) we have the radial distribution function which for the hydrogen ground state is obtained by multiplying the square of the wavefunction by a spherical shell volume element (\(4\pi r^2\)).

\[ \text{Radial distribution function}(r) = |\Psi(r)|^2 \cdot 4\pi r^2 \]

We will use atomic units \(a_{0}=1,\hbar=1,me=1,e=1\)

Notice that because of the units the probability is not normalized (the values do not integrate to one)

#Lets plot the radial distribution functions

r = np.linspace(0,10,100)
n=1
l=0
R = psi_R(r,n,l)
plt.plot(r, 4*np.pi*(r**2 )* R**2, lw=3)
plt.xlabel('$r [a_0]$',fontsize=20)
plt.ylabel('$4\pi r^{2} \cdot R^2_{nl}(r)$', fontsize=20)
plt.title("Radial probability distribution for n="+str(n)+" l="+str(l))
plt.grid('True')

#Comparing radial distribution functions
fig = plt.figure(figsize=(12,7))
ax1 = plt.subplot(121) #121 means: the subplots will be 1 row 2 columns and this one is subplot number 1
ax2 = plt.subplot(122, sharey = ax1) #share-y means that it will share the y-axis with the ax1 plot

r = np.linspace(0,10,100)
n=2
l=0
R = psi_R(r,n,l)
ax1.plot(r, 4*np.pi*(r**2 )* R**2, lw=3)
ax1.set_xlabel('$r [a_0]$',fontsize=20)
ax1.set_ylabel('$4\pi r^{2} \cdot R^2_{nl}(r)$', fontsize=20)
ax1.set_title("Radial probability distribution for n="+str(n)+" l="+str(l))
ax1.grid('True')
n=2
l=1
R = psi_R(r,n,l)
ax2.plot(r, 4*np.pi*(r**2 )* R**2, lw=3)
ax2.set_xlabel('$r [a_0]$',fontsize=20)
ax2.set_ylabel('$4\pi r^{2} \cdot R^2_{nl}(r)$', fontsize=20)
ax2.set_title("Radial probability distribution for n="+str(n)+" l="+str(l))
ax2.grid('True')

or more interactively

nmax=10
import numpy as np
import ipywidgets as widgets
@widgets.interact(n = np.arange(1,nmax,1), l = np.arange(0,nmax-1,1))

def plot_radial(n=1,l=0):
    r =    np.linspace(0,250,10000)
    psi2 = psi_R(r,n,l)**2 * (r**2) * 4*np.pi
    
    plt.plot(r, psi2, lw=2, color='red')
    plt.xlabel('$r [a_0]$')
    plt.ylabel('$4\pi r^2 R^2_{nl}(r)$')
    rmax = n**2*(1+0.5*(1-l*(l+1)/n**2))
    plt.xlim([0, 2*rmax])

• What is the distance that an electron in a 1s orbital most probable to be? And in 2s? And in 2p?

• Are the most probable distance the average distance? The average distance is called expectation value. Would you expect the expectation value of the 1s be higher or lower than the most probable distance?

• Predict if the expected value for “r” of 2s is lower or higher than for 2p.

Plotting the entire wavefunction: probability and cutoff

If we want to plot the entire wavefunction we need to plot that huge function that depends on \(r\), \(\theta\) and \(\phi\). Luckily for us, again, all the polynomials (the Spherical harmonics and Laguerre) are available on python. Even though the angular function depends on \(\theta\) and \(\phi\) it will be more visually understandable if we plot the wavefunctions as a function of x,y,z So rather than plotting it as a function of the angles, we will build a mesh of xyz and calculate the value of the wavefunction at that point.

One problem is that this is a 3D function, so we would need a fourth dimension to plot its value. The only solution is to plot an “isosurface”, that is, plot a surface not a volume where the function has a specific value. This is why we need to come up with a value. The variable iso is set as 30% of the maximum value of the function.

First we will calculate the whole wavefunction for a given n,l, and m.

#source: https://github.com/wmfschneider/CHE30324/blob/master/Outline/CHE30324-outline.org
import plotly.graph_objects as go
from sympy.physics.hydrogen import R_nl
from sympy import Ynm
from sympy import var

var("r theta phi")

n = 3
l = 2
m = 0
Z = 1 # this is the charge of the nucleus. Make it Z=2 to see He+ or Z=3 for Li(2+)

xmesh=[]
ymesh=[]
zmesh=[]
psi=[]
points=20

for z in np.linspace(-10.,10.,points):
    for y in np.linspace(-10.,10.,points):
        for x in np.linspace(-10.,10.,points):
            xmesh.append(x)
            ymesh.append(y)
            zmesh.append(z)
            r = np.sqrt(x**2+y**2+z**2)
            theta=np.arccos(z/r)
            phi=np.arctan(y/x)
            value=R_nl(n,l,r,Z)* Ynm(l,m,theta,phi).evalf() 
            psi.append(float(value))

and now we can plot it deciding on the value of iso

#source: https://github.com/wmfschneider/CHE30324/blob/master/Outline/CHE30324-outline.org
iso = 0.2*max(psi)

fig= go.Figure(data=go.Isosurface(
    x=xmesh,
    y=ymesh,
    z=zmesh,
    value=psi,
    isomin=-iso,
    isomax=iso,
    surface_count=2,
    opacity=1.0,
    caps=dict(x_show=False, y_show=False,z_show=False)
))

fig.show()

Quantum weirdness

Popular (non-academic) Resources:

• On Schrodinger’s cat and quantum weirdness: https://www.youtube.com/watch?v=uWMTOrux0LM

• The double slit experiment: Quantum superposition, the wave nature of reality and the role of observation https://www.youtube.com/watch?v=NvzSLByrw4Q

Critical thinking questions on quantum weirdness

It seems that “reality” exists as many situations at the same time. In other words, it is not that the electron is a wave, rather, it is that we can only know about quantum objects as a wave (superposition) of possibilities, and it is not until we measure the system that one of those possibilities is defined. In other words, our observing changes the nature of what we observe. Can objective science exist? If the observer cannot observe a system without changing it, is there an independent reality outside of the observer?

More philosophical implications here: https://plato.stanford.edu/entries/qt-issues/

Exercise

• Plot the energy levels for the atom of hydrogen for each allowed combination of n,l, and m.

• What is the degeneracy level for each n?

• Plot the radial wavefunction for the 2p and 3d.

• Plot the radial distribution functions for 3s and 3d. What is the most probable distance to the nucleus for these two orbitals? Do you think the most probable is the average distance?

• Is the most probable distance the size of the atom?

• Use plotly to plot the entire wavefunction for 2 orbitals of the 3d. What cutoff or iso value did you use? What’s a cutoff and why do we need it?